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Analysis of slow transitions between nonequilibrium steady states 
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Transitions between nonequilibrium steady states obey a generalized Clausius inequality, which 
becomes an equality in the quasistatic limit. For slow but finite transitions, we show that the 
behavior of the system is described by a response matrix whose elements are given by a far-from- 
equilibrium Green-Kubo formula, involving the decay of correlations evaluated in the nonequilibrium 
steady state. This result leads to a fluctuation-dissipation relation between the mean and variance 
of the nonadiabatic entropy production, As na . Furthermore, our results extend - to nonequilibrium 
steady states - the thermodynamic metric structure introduced by Sivak and Crooks for analyzing 
minimal-dissipation protocols for transitions between equilibrium states. 


Classical thermodynamics is a macroscopic theory 
built around the concept of the equilibrium state [l|, 
whereas statistical mechanics is a microscopic theory that 
represents equilibrium states in terms of known statisti¬ 
cal ensembles Q- Both frameworks accurately describe 
the equilibrium properties of matter and the constraints 
that must be satisfied when systems undergo transitions 
between equilibrium states. No comparably general the¬ 
ories exist for systems away from equilibrium, not even 
for steady states. We have no a priori microscopic rep¬ 
resentations of nonequilibrium steady states; neither do 
we have a macroscopic, empirical thermodynamic theory 
to compare to. 

A phenomenological thermodynamic theory of 
nonequilibrium steady states was suggested by Oono 
and Paniconi [sj, who introduced the crucial concept of 
housekeeping heat , which is the heat transfer necessary 
to maintain a system in a given steady state. The excess 
heat Q ex is then the renormalization of the total heat Q 
by the housekeeping heat QhG 

Qex — Q Q hk- (1) 

In the context of Markovian dynamics, these notions were 
modified somewhat and given quantitative expressions by 
Hatano and Sasa 0, Speck and Seifert jU and Ge and 
Qian B0, leading to derivations of a generalized Clau¬ 
sius inequality for transitions between steady states [d, bj: 

AS +J dtp(t)(Qe X ) > o. (2) 

Here AS denotes the net change in the entropy of the 
system, and (Qex) is the ensemble-averaged rate of excess 
heat transfer to a reservoir at inverse temperature /3(f). 
(See also Bertini et al Si for a related inequality for 
renormalized work.) Under equilibrium dynamics Eq. [2] 
reduces to the usual Clausius inequality, as Qhk = 0 and 
hence (f @x — Q • 

Eq. [2] suggests that there may exist a steady-state ther¬ 
modynamic framework that closely parallels equilibrium 
thermodynamics, with steady states and excess heat in 


the roles traditionally assigned to equilibrium states and 
heat, respectively. In this Letter we further develop these 
parallels, by analyzing finitely slow transitions from one 
nonequilibrium steady state to another. At leading order 
in a perturbative expansion in the driving speed, Eq. [2] 
becomes an equality and hysteresis vanishes, suggesting 
that quasistatic steady-state processes are natural coun¬ 
terparts of reversible thermodynamic processes @]. At 
the next order, a response matrix £ governs the renormal¬ 
ized entropy production (Eq. [20]). We show that the ele¬ 
ments Qjy are given by time-integrated correlation func¬ 
tions evaluated in the steady state (Eq. l23|) . in exact 
analogy with equilibrium Green-Kubo relations. These 
results lead to a fluctuation-dissipation theorem, which 
states that the average renormalized entropy production 
is equal to half its variance (Eo. l25l) . Finally, our analy¬ 
sis generalizes recent progress related to thermodynamic 
length and its applications to the determination of opti¬ 
mal driving protocols loj- 13 • 

Our results reveal that the stochastic theory of slow 
transitions between equilibrium states [Tjl El, 151 ex¬ 


tends directly to slow transitions between nonequilibrium 
steady states, under the definition of housekeeping heat 
proposed in Refs. H0- We note that alternative defini¬ 
tions have been suggested, and corresponding generalized 
Clausius inequalities have been derived, by Komatsu et 
al 16, 171 and Maes and Netocny [Til ]. It remains an open 
question whether the results we derive have counterparts 
within the frameworks of Refs. EMI- 

We develop our theory within the context of isothermal 
Markovian dynamics on a network, in which the system’s 
evolution consists of random transitions among a set of 
discrete states. This mesoscopic level of description is 
well-suited for small stochastic systems such as molec¬ 
ular motors 0. We expect that our results extend as 
well to diffusive processes, described by Langevin [20j ] 
and Fokker-Planck equations [21]. In the following three 
paragraphs we specify notation and define quantities that 
will play important roles in our subsequent analysis. 

Consider a system with N discrete states i £ 
{1,2,... N}, in contact with a reservoir at temperature 
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/3 _1 . Induced by thermal fluctuations from the reservoir, 
the system makes random, Poissonian transitions among 
its states, with R^ > 0 denoting the transition rate from 
j to i. The probability distribution p = (pi,p 2) • ■ -Pn) T 
(the superscript T denotes transposition) satisfies 

|p = Kp, (3) 

where the rate matrix TZ is formed by the rates Rij. The 
diagonal elements satisfy Ra = — , • Rji to preserve 

normalization: yu Pi(t) = 1. The current 

Jij Fij Pj RjiPi (4) 


is the instantaneous flow of probability from j to i. 
The quantity Qij = /3 _1 In (Rij/Rji) represents the heat 
transferred from the system to its thermal surroundings, 
during a transition from j to i [22| . Thus the average 
instantaneous rate of heat transfer to the reservoir is 


(Q) = r 1 E J b ln (§ £ )- 

i<j ^ / 


(5) 


Consider now a set of parameters A = {Ai, A 2 ,... Xk } 
that determine the transition rates R%j{ A) and corre¬ 
sponding steady states t r A . The generalized Clausius 
inequality, Eq. [21 applies to processes in which the sys¬ 
tem is driven from state 7r A to state 7t b , by varying 
these parameters from A to B over a time interval At. 
We assume that the protocol A(t) is smooth, so that 
d\/dt is well-defined. The term AS appearing in Eq. [2] 
is the change in Shannon entropy, S^b) — S( tv a ), with 
S = i 7q In 77^. 

We will analyze Eq. [2] for processes in which the pa¬ 
rameters are varied slowly , hence the system remains near 
the nonequilibrium steady state. For convenience, we in¬ 
troduce a small parameter e oc \d\/dt\, so that At oc e -1 
for fixed A and B. We will find that in the quasistatic 
limit (e —> 0) as well as in the leading correction to this 
limit, the response of the system is remarkably analogous 
to that of a slowly driven equilibrium system. 

Let us write p(t) as the sum of the instantaneous 
steady state distribution Tv(t) = 7r A( ^ and a small cor¬ 
rection, or “lag”, Ap(f): 


We further assume that (i) Rij ^ 0 if and only if Rji ^ 
0, and (ii) any state i can be reached from any other 
state j either directly or via intermediate states. Under 
these assumptions, any initial distribution p(0) relaxes 
to a unique steady state distribution 7T, with TZtv = 0, 
characterized by steady currents J?? = RijTVj—RjiTVi [2d| ■ 
If Jfj = 0 for all i j, then the dynamics generated by 1Z 
satisfy detailed balance, and the distribution 7r represents 
an equilibrium state. If some ^ 0 - as we generically 
assume throughout this Letter - then detailed balance is 
broken and tv specifies a nonequilibrium steady state. 

We now define the quantities 


Fij = In 


f Fij n j \ 

V Rjiivi J ’ 


( 6 ) 


which we interpret as effective thermodynamic forces 
that drive the probability currents in nonequilibrium 
steady states (note that F?? = 0 if and only if = 0). 
These forces are uniquely determined by the rate matrix 
7 Z, and are nonlocal, in the sense that each F?? gener¬ 
ally depends on all the elements of 7 Z, via the stationary 
distribution 7r. The instantaneous rate of housekeeping 
heat transfer to the medium is now defined as \j] 


(Qhk)=/3- 1 ^J^, (7) 

i<j 

which has a natural interpretation as the power associ¬ 
ated with the forces F??. It can be shown that (Qhk) > 0 
for any distribution p, and ( Q) = (Q h k) in the steady 
state p = 7T 7]. Using Eq. [I] we obtain the excess 
heat transfer rate, 


p(7) = tv (t) + Ap(t). 


(9) 


Substituting Eq. [9] into Eq. [3] and using the condition 
TZtv = 0 , we obtain the linear inhomogeneous equation, 


d d 

-Ap-KAp = --*. 


( 10 ) 


Let us_ now define a generalized inverse TZ + by the rela¬ 
tions 


23, 


(TZTZ + ) ij - (TZ + TZ) ij = 5 ij -Tv i (lla) 

1Z+TV = 0 , 1 t TZ + = 0 T , (lib) 

where Sij is the Kronecker delta, 0 T = (0, 0, - • ■ 0), and 
l r = (1,1, ■ • -1). Applying TZ + to both sides of Eq. [lOl 
we get 


1 - 7W 


dt 


ap = K+ 5U 


which is solved iteratively to obtain 

d 

n—1 


Ap = J2\K 


dt 




3n£ 


( 12 ) 


(13) 


n—1 


This expansion of the lag Ap in powers of the driving rate 
leads to a corresponding expansion ( Q ex ) = b n e n 

(via Eq. 1141 belowL In the following, we examine in turn 
the first two terms of this expansion. 

We start by rewriting Eq. [S] in the form [f| 0] 


(Qex) = 1 E R SPi ln7r i- ( 14 ) 

hj 

From Eq. [9] and the n = 1 term of Eq. [13] we get 


(Qex) = (Q) - (Qhk) = P 1 E Ji 0 111 

i<j 


( 8 ) 


(l (Qex) — ^ ( R j j TV j 111 TV, 


^2 Rij Rf k TVk In 1 


(15) 


i,j,k 
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The first term on the right vanishes, since TZtv = 0. Us¬ 
ing Eq. Illal in the second term, we get a sum of two 
terms: (1) JA 7Tj In 7r ? ;, which is equal to — dS/dt , and 
(2) — k TTiTTk In7 Ti, which vanishes by conservation of 
normalization: 7f*, = 0. We thus arrive at 

m ex ) = -§ + G(e% (16) 

which implies that in the quasistatic limit (e —» 0, with 
At oc e -1 ), Eq. [2] becomes an equality: 


obtain, after some simplifying steps (see SI), 


A S + J d t/3(Q ex ) 

= -A V In- - 




i,3 


i,3 


dt 13 dt 


(19) 


If we now assume that dX/dt = 0 at the start and end 
of the process, then the first term on the right of Eq. fl9l 
vanishes. As the steady states n are determined by the 
parameters A, we can rewrite Eq. M in the form 


AS +J dtm(Qe X ) = 0 (17) 


This result is a generalized Clausius equality for qua¬ 
sistatic transitions between nonequilibrium steady states. 
An equivalent result was obtained for overdamped 
Langevin processes in Ref. 0,|I3|- Eci. H71 implies that the 
integral f dt f3(t)(Q ex ) is independent of the quasistatic 
path taken from A to B in A-space, and therefore van¬ 
ishes when the path is cyclic. (Interestingly, if Qhk is 
defined as in Refs. 3, 17], then for cyclic paths this in¬ 
tegral is described in terms of a geometric phase [25| .) 


Recall that reversible equilibrium processes, which sat¬ 
isfy AS + f d t/3(Q) = 0, are characterized by zero en¬ 
tropy production in the universe: any change in the sys¬ 
tem’s entropy is balanced by a compensating change in its 
surroundings. By analogy, in quasistatic nonequilibrium 
processes, which satisfy Eq.El the entropy change of the 
system, AS , is balanced by the excess entropy produced 
in the reservoir, j d t/3(Q ex ). (The total entropy pro¬ 
duction in the reservoir diverges in the quasistatic limit, 
f dt/3(Q) —> oo, due to the continual flow of housekeep¬ 
ing heat.) Moreover, just as a system remains arbitrarily 
close to equilibrium during a reversible processes, a sys¬ 
tem undergoing a quasistatic nonequilibrium transition 
remains arbitrarily close to the nonequilibrium steady 
state (Ap oc e). In both cases, equilibrium and nonequi¬ 
librium, the system retraces its path in the reverse order 
when it is subjected to the reverse process A : A «— B; 
in this sense, there is no hysteresis. In view of these par¬ 
allels, it is natural to think of quasistatic nonequilibrium 
processes as the nonequilibrium analogues of reversible 
equilibrium processes, as suggested by Oono and Pani- 
coni [3j. 


Let us now move beyond the quasistatic limit, by in¬ 
cluding the n = 2 term of Eq.[l3]in the analysis. Starting 
with Eq. 0 we obtain 


/3(Qex) = ( 18 ) 

i,j k,l 


in place of Eq. 1161 Integrating with respect to time, we 


A S + j df/3(Q e x} =Jdt\ £(A)A 

= / dt a t 'c(A)a, 


( 20 ) 


where ( = (£ + £ T )/2 is the symmetric part of a matrix 
£(A) whose elements are 


V- 9 In 7 u , d In 7r, 

~ > 7r,-—- RZ- 




d\„ 


13 d\„ 


( 21 ) 


Equation [2D] provides the leading correction to Eq. El 
and is the counterpart of analogous results for slow tran¬ 
sitions between equilibrium states [3 EM!. 

We now derive a Green-Kubo relation for the elements 
of the matrix £(A). Let us define a set of observables 


FtW = 


91n7Tj(A) 

dX u 


d = 1, 


• K. 


( 22 ) 


When the system is in the steady state 7t a , its microstate 
i(t) € {1,2,... N} fluctuates in time, hence so does each 
F M (t) = Ffi t y around a mean value (F ll )\ = 0. Letting 
(F tJ, (0)F v (t))\ denote a correlation function evaluated in 
the nonequilibrium steady state, the matrix elements 
can be rewritten as (see SI for details): 


Cm-(A) - 


1 


/*+oo 


dt (F li (0)F v {t)) 


A- 


(23) 


This result relates an excess dissipation coefficient Cm- 
to stationary fluctuations in the nonequilibrium steady 
state. (Analogously, for near-equilibrium transitions 
the friction tensor is determined by equilibrium fluctua¬ 
tions 0 .) We emphasize that the steady state in Ea.l23l 
may be far from thermal equilibrium. 

As shown by Prost et al [3|, and for general Markov 
processes by Hanggi and Thomas [13] , an expression sim¬ 
ilar to Eq. [23] describes the linear response of a system 
to small perturbations around a given steady state. By 
contrast, our analysis applies to slow transitions between 
two steady states that may differ substantially. 

The left side of Eq. [2] (or Eq. I2P1) is the ensemble aver¬ 
age of a quantity identified by Esposito et al 2^, 2 3 as 
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the nonadiabatic component of entropy production: 

As na = - In 7 t® t) + In 7r^ 0) + J dt f3Q ex , (24) 

where i(t) is the microstate of the system during a single 
realization of the process. For slow driving the mean and 
variance of As na satisfy (see SI) 

(As na ) = —a A Sna - (25) 


of the Fisher information matrix [32| of the steady state 
distribution tt x , 


V( A ) = 5>i (5 M ln7rO (d v hun) = (W) A , (30) 


and a matrix generalization r(A) of the relaxation time, 


_ r,,wi o)k 

^ J 0 (F»F") X 


(31) 


Let us place this result in context. If the dynamics sat¬ 
isfy detailed balance and temperature is constant, Eg. HH1 
reduces to As na = (3{W — A F), where W is the work 
performed on the system and A F is the free energy dif¬ 
ference between equilibrium states A and B. Ea.l25lthen 
becomes 


{W)=AF+^a 2 w . 


(26) 


This fluctuation-dissipation relation for isothermal, near¬ 
equilibrium processes was originally proposed by Her¬ 
mans (3Q) and Wood 311, using a Gaussian assumption 
for the work distribution, and more recently has been ob¬ 
tained using projection operator techniques by Speck and 
Seifert 0 and multiple time-scale analysis by Hoppenau 
and Engel [15:]. Our result, Eq. [25] generalizes Ea. l26lto 
slow transitions between nonequilibrium steady states. 
Using the techniques of Refs. 14,[l5|, one can show that 
for such transitions, typical values of As na follow a Gaus¬ 
sian distribution. 

Since the matrix £(A) is positive semidefinite (see SI), 
it provides a metric in A-space, related to nonadiabatic 
entropy production for slow processes. Specifically, for a 
protocol A(f) describing a contour C in parameter space, 
we define a length lc by the contour integral 

lc = l \J dA T £(A) dA. (27) 

By the Cauchy-Schwarz inequality we have 


dtA T C(A)A> J^, 


(28) 


where At is the duration of the transition. If the driv¬ 
ing is sufficiently slow that Eq. [221 applies, then the right 
side of Eq. EBlprovides a lower bound for the nonadiabatic 
entropy production: (As na ) > l^/At. This result gener¬ 
alizes a similar structure for transitions between equilib¬ 
rium states, which has been used to determine optimal 
(minimally dissipative) protocols A(f) opt for fixed end¬ 
points A and B and duration At 10Ml3l|. 

Finally, following Sivak and Crooks [12], we note that 
the metric £(A) is the Hadamard (term by term) product 


(29) 


In the framework of finite-time thermodynamics 
built on the metric geometry of equilibrium thermody¬ 


namics [36 4l|, the Fisher-information matrix X(A) alone 
is assumed to dictate the average dissipation of heat 
during a thermodynamic process. The current work, in 
agreement with Refs. [12|, [13j, [42], emphasizes the impor¬ 
tance of relaxation dynamics, via the matrix r(A). 

This material is based upon work supported by the 
U.S. Army Research Laboratory and the U.S. Army Re¬ 
search Office under contract number W911NF-13-1-0390. 
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SUPPLEMENTARY INFORMATION 


We begin by deriving Eq. IS 101 below, which we will 
later use to derive Eq. 23 of the main text. We note that 
Eq. 11, which uniquely specifies the matrix 7?. + , 1 can be 
written as a set of four equations: 


1ZTZ + 

= I-tv1 t 

(SI) 

TZ+TZ 

= I — tv1 t 

(S2) 

n + tv 

= 0 

(S3) 

1 T 1Z+ 

= 0 T 

(S4) 


where I is the identity matrix. In the following, we show 
that the integral 

pOO 

M= dte nt (Tvl T -I) (S5) 

Jo 

satisfies all these defining relations of 7Z + , thereby con¬ 
stituting an exact expression of the latter. In these cal¬ 
culations, 7Z and 7r are fixed (i.e. time-independent). 

Eq. E2 Multiplying M by 7Z on the left, we get 



The third line follows from the relaxation property: 
lim^oo exp (TZt)p(O) = 7r for any normalized p(0), hence 
lim^oo exp (7 Zt) = 7rl T . The last line follows from the 
normalization l 1 tv = 1 and a cancellation of terms in the 
product. 

Eq. El Multiplying M by 1Z on the right, we get 


M1Z= / dte'^inV -1)7?. 


pOO 

/ dt e nt lZ(Tvl T — I) 

J 0 

p£=00 

/ d(e™)(7rl T -I) 

Jt =0 


(S7) 


H—0 

= I — tvI 1 . 


The second line utilizes the relations: l r 7Z 
TZtv = 0. The last line follows from Eq. IS6I 
Eq. E3 Multiplying M by 7r we get 

POO 

A4tv = / dt e nt (ivl r — I)7r 

Jo 

dt e nt (iv — 7r) = 0. 



0 T and 


(S8) 


The second line follows from the normalization condition 

1 T 7T = 1. 

Eq. Multiplying M. by 1 T we get 

POO 

1 t M= / dtl r e TCt (7rl T -I) 

Jo 

dtl T (7Tl T ^I) (S9) 

dt(l T - 1 T ) = 0 T . 

The second line follows from the stationary relation 
l T e 727 = l 7 . The third line utilizes normalization condi¬ 
tion of tv. Comparing Eqs.[S6llSl]with Eas. fsTlFvtl we see 
that Ai satifies all the defining relations of 7? + , leading 
to the identification 

POO 

1Z + =M= dte nt (Tvl T — I). (S10) 

Jo 



Derivation of Eq. 19 


We start with Eq. 18: 

/3(Qex> = ~ + E R d ( R tl*l) ■ ( S11 ) 

ijkl 


Summing over the index j and using Eq. 11a, the second 
term on the right becomes 




- 7Ti) 





(S12) 


Summing over i and using S = — JT tt,; In 7q, Eq. IS12I 
becomes 

B^TTfc +S)4: ■ ( S13 ) 


The second term in Ea. IS 131 (containing S) vanishes upon 
summing over k, since 1 T 7Z + = 0 T (Eq. lib). Combining 
results, we get 

/3(<0ex) = + B ln7rfe ^ ■ ( S14 ) 


Integrating both sides with respect to time leads us to 
Eq. 19, after integration by parts. 


Derivation of Eq. 23 

From the definition of (Eq. 22), we have 


1 In fact, there is some redundancy: given Eq. 11a, the two equa¬ 
tions appearing in Eq. lib are equivalent to one another. 
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with all quantities evaluated at fixed A. We can rewrite 
the right hand side as 
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where in the third line we have used l 7 ■ (dir/dXP) = 0, 
which follows from the normalization condition l r 7r = 1; 
in the fourth line we have used the stationarity condition 
e nt ir = 7r; and in the last line we have used Eq. IS10I 
Combining Eos. IS15l and lS16l we get 
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using Eq. 21. We also have 
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(S18) 

using Ea. lS17l Combining Eas. IS17I and IS18I gives us the 
desired result: 

U = \ + U) = \ £°dt (F^0)F"(t)) x . (S19) 


Derivation of Eq. 25 

By Eq. 8, the excess heat associated with a transition 
from state j to state i is equal to /3 -1 \n{in/irj). Thus the 
total excess heat dissipated during a single realization of 
the process is given by 
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where i(t n —) and i[t n +) denote the states of the system 
just before and just after a transition at time t n . Eq. 24 


now becomes, with t as the duration of the process, 
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where F^(t) = Ftt^ (A(f)). In going from the first line 
to the second, we have used the fact that the net change 
in the value of hi7r^y \ from t = 0 to t = r, is a sum 
of contributions due to: (1) discrete changes in i(t) at 
the transitions, and (2) the continuous variation of A(t) 
between transitions. 

Over an ensemble of realizations, we have 


((As na ) 2 ) = dt df y2x^t)X^t')(F^t)F^t')). 
do do 

(S22) 

Let us now evaluate this expression in the limit of slow 
driving, where there is a separation of time scales be¬ 
tween the slow variation of A(f) and the fast evolution 
of F(i) (due to rapid transitions between states). As in 
the main text, it is convenient to think in terms of a 
small parameter e oc |A|, so that changes in A(f) occur 
on time scales of order e -1 , while changes in F(t) occur 
on time scales of order unity. The correlation function 
(F ^t (t)F ^ '(f , )) decays with \t' — 1\, on a time scale of order 
unity. As a result, to leading order we can replace this 
correlation function with one that is evaluated in a given 
nonequilibrium steady state: 


{F^{t)F v {t')) » {F»{t)F v (t')) X = C^(s; A), (S23) 


with s = t' — t and A = A(t) in Eq. IS23I Here as in 
the main text, F^(t) is evaluated along a trajectory gen¬ 
erated at fixed A, unlike - defined above - which 

is evaluated along a trajectory evolving under the slow 
variation of the external parameters. Eq. IS22I now be¬ 
comes 
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Since r oc e -1 , we conclude that for most values of t 
between 0 and r, both t and T — t are much larger than 
the time scale over which the correlation function decays. 
Hence to leading order in e we can write 
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using Eq. IS19I Now note that Eq. 20 can be written as 

<As na ) = [dt\ T a\)\ . (S25) 

Jo 

Since A oc e and r oc e” 1 , both ((As na ) 2 ) and (As na ) 
scale as e, which implies that 

°As na = ((A S „ a ) 2 } - (As na ) 2 » ((As na ) 2 ). (S26) 

Hence to leading order we have (As na ) = (1/2 )ct^ . 


Squaring the value of Y and averaging over an ensemble 
of steady-state trajectories, we get 

pAt pAt 

(Y 2 ) = / dt / At'(F^(t)F v (t'))x 

pAt pAt-t 

= / At / dsC M!/ (s;A) 

/*+A t 

= V a M a, / ds (At - |s|) C^(s; A) (S28) 


If we now divide both sides by At and consider the limit 
At —> oo, we arrive at 


The matrix £(A) is positive semidefinite 

To establish this result, let i{t) denote a trajectory 
evolving under the stationary dynamics at fixed A, let 
At > 0 be an interval of time, let (ai, 02 , ■ • ■ a k} denote 
a set of real values, and consider a quantity 

pAt 
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Since the left side is non-negative for any choice of 
{ai, 02 , • • • a^}, we conclude that £(A) must be positive 
semidefinite. 





